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Abstract 

We present a stable and convergent method for solving initial value problems based on 
the use of differentiation matrices obtained by Lagrange interpolation. This implicit 
multistep-like method is easy-to-use and performs pretty well in the solution of mildly 
stiff problems and it can also be applied directly to differential problems in the complex 
plane. 
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1 Introduction 



A general technique for obtaining approximations for the derivative of a well-behaved 
function x(t) consists in expressing the derivative x'(t) evaluated at an arbitrary point 
as a linear combination of the function values x(tk) at iV + 1 nodes tj., k = 0, 1, . . . , N, 
i.e., x'(t) ~ Yl!k=Q^k x ^k) wnere the numbers 7^ must be such that the truncation error 
tends to zero as the mesh goes to zero [1] . This N + 1-point differentiation scheme can 
be realized by using the Lagrange interpolating polynomial. In spite of the fact that 
Lagrange interpolation has been widely and greatly used in numerical differentiation, it 
seems that a very simple and powerful simplification has been overlooked for many years. 
If the N + 1-point derivative of x(t) obtained by Lagrange interpolation is evaluated at 
the nodes, it can be written the form [21 [3l H] 



N 

AW = E^**) + 7ivTTv p/(tj)x(7V+1)(rj) ' (1) 

fc=0 ^ >' 

N 



where Tj G (t Q ,t N ), P(t) = H k=0 {t ~~ and 



T> ik = < 



( N 1 



p'it 



k [tj — tk) P'(tk) 



(2) 



If the function x(t) to be differentiated is a polynomial of degree at most N, the truncation 
error is zero and yield the exact values x'(tj) from which the function x'(t) can be 
retrieved by an interpolation. This is why the differentiation matrix D is a projection 
of d/dt in the subspace of polynomials of degree at most N. The arbitrariness of the 
nodes and the matrix form of Eq. (jTJ), x' = T>x + e, where T> is the matrix whose 
elements are given by (j2J) and e is the vector corresponding to the truncation error, 
have given rise to a simple method for finding accurate numerical solutions to two-point 
boundary value problems [5] with a relatively small number of nodes. However, the 
accuracy achieved by this technique depends on the selection of an auxiliary function 
needed to incorporate the boundary conditions. On the other hand, the definition (j2J) 
can be generalized straightforwardly to give a differentiation matrix for meromorphic 
functions, yielding a method for solving singular differential problems in the complex 
plane [1]. 

Our aim in this paper is to implement a simple method to solve the initial-value problem 
by using the differentiation matrix ([2]). This is done in the next section. It is shown that 
this method is convergent and stable and some numerical tests with benchmark problems 
are given. 
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2 Outline of the method 



Consider the IVP 

x'(t) = f(x, t), x(a) = a, (3) 

defined on 71 = {(t, x)\ a < t < b, c < x < d}, where f{x,t) is Lipschitz on TZ in the 
variable x. Let ai be a point of (a, b) and evaluate the differential equation of ([3]) at the 
N + 1 nodes a = t < t\ < t 2 ■ ■ ■ < t N = a±. This yields the vector equation 



/ x'{a) \ 
x'(t 2 ) 



a, a 



/ /( 

f(X2,t 2 ) 



\ 



\x'(t N )J \f(x N ,t N )J 
where Xk stands for x(tk). According to ([!]), this equation can be approximated by 

/ f(a, a) \ 
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(4) 



This system has only N unknowns £i, . . . , £tv> that can be found by solving 

N 



J2 V i& - f&> *i) = ~ ad j> ^ = 1, 2, • • • , N, 



(5) 



k=l 



where dj = Vj a . This yields approximations to x(tj) in [a, a\]. Now, let us define 
«i = £,n and choose a 2 G (a\, b]. Then, the local problem 

y'(t) = f(y,t), te[a 1 ,a 2 ], y(a 1 )=a 1 , 

can be solved numerically along the above lines to give new approximations £j to x(tj) 
in [01,02]. Defining a2 as the new £at, this procedure is repeated from the subinterval 
[a n _i, a n ] to the subinterval [a n , a n+ x] until the final point b is reached. 



3 Main differences with popular IVP methods 

It must be noted that the proposed method, from now on referred to as the Simple Conver- 
gent Solver (SCS), is based on the system of equations (Q. Since the solution of the latter 
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produces the approximations to the unknown values x(t\), x(t 2 ), x(t]^) simultaneously, 
SCS is an implicit method. It looks like a linear multistep method, however, the fact 
that N unknowns are calculated simultaneously in every step makes SCS different from 
those kind of methods. On the other hand, even though SCS is based on Lagrange inter- 
polation and incorporates N values f(£j,tj), the use of differentiation matrices, instead 
of quadratures, is what makes SCS also different from the general multistep-multistage 
methods discussed in the literature [t^ [TUj. 



4 Analytic properties 

Note that the system (jSJ) has the matrix form 

Dnin - = -a n -id n , (6) 

where the index n indicates that the variables are being considered in the subinterval 
[a n _i, a n ]. In order to simplify the notation, let us rewrite this equation as 

D£-ft= -ad, (7) 

where = £ nj , (f^j = f(£ n j, t n j), dj = d nj = V ja , and D = (D jk ) is the N x N submatrix 
of V with Dj k = Vjk, j, k — 1, 2, . . . , N, and Vj k being computed according to (|2J) at the 
nodes a n _i = t n0 < t n \ < t n2 < . . . < t n ^ = a n . 

As mentioned above, the approximation to the solution of the IVP in [a n _i,a n ] is 
obtained by solving the system of equations (jTj). In general, for a nonlinear f(x,t), a 
Newtonian iteration can be applied in order to solve the equations. It must be noted, 
however, that this is not a major problem: the required Jacobian can be easily approx- 
imated, if needed, by using standard differences due to the simple form of the left-hand 
side, since only the partial derivative df(x, t)/dx is required. Thus, the (k + l)th iteration 
can be written as 

where A^ fe ^ is the diagonal matrix whose non-zero elements are 

fdf(^,t nl )/dx\ 

df{^ 2 \t n2 )/dx 

\df{^lt nN )/dx) 

and is the vector with elements (f^)j = f(£,nj,t n j)- The initial trial can be 
taken as the vector with all its components equal to the initial condition at a n _i. 
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On the other hand, if / is linear, let us say f(x,t) = kx + (j)(t), the approximation £ is 
the solution of (j7j) which becomes 

(D — kIn)£, — —ad — <p, 

where In is the identity matrix of dimension N and is the vector of entries <j)(t n j), 

.1 l,2 v. 

It is clear that the existence of the numerical approximation given by the method is based 
on the invertibility of the sum of D and a diagonal matrix. 

The invertibility of D can be easily proved. Let be the differentiation matrix con- 
structed according fl2]) with the set of N points t n \ < t n2 < ■ ■ ■ < t nN (note that the first 
point t n0 has been removed). Then a little algebra shows that 

(T - t n0 l N )D = (T — t n0 l N )V N + In, (9) 

where T is the diagonal matrix whose non-zero entries are t n \, trcii • • • , t n j\f. Since T^n is a 
differentiation matrix, it yields exact values for the derivatives of a polynomial of degree 
m, m < N — 1 [2J. Therefore, the derivatives of (t — t n o) m , m — 0, 1, . . . , N — 1, can be 
reproduced for each n applying to the vector of entries (t n j — t n0 ) m , j = 1, 2, . . . , N. 
Thus, we have that 
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The substitution of this result in (Q yields 
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This equation shows that the eigenvalues of (T — t n Ql^)D are the integers 1, 2, . . . , N for 
any set of points t n \ < t n i < ■ ■ ■ < t n N- Therefore, (T — t n0 i7v)-D is invertible. Since 
(T — t„oljv) is invertible, we have proved the 

Lemma 4.1. D is invertible. 

This lemma will be of great importance to prove convergence and stability for the 
method. More properties of the matrix D will be studied in more detail elsewhere. 
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5 Convergence, consistency and stability 



As noted above, the proposed method is different from the standard methods for initial 
value problems, since it works as a "block-implicit" method. Thus, it is convenient to 
discuss ad hoc proofs of its basic properties. 
Let us assume that the solution y(t) of the local problem 

y'{t) = f{y,t), £e[a n _i,a„], y(a n _ x ) = a n _i, 

is sufficiently smooth and Lipschitz in y. 

To prove convergence, let us denote by y = (y n i, y n 2, • • ■ , UnN) T the vector whose entries 
are the values of y(t) at the N nodes t n i, t n 2, . . . , t n u, i.e., y n j = y(t n j), and assume that 
[a, b] has been divided in M subintervals. Then, applying equation ([1]), we get 

N E ■ 

^ DjkVnk ~ f(Vnj, t nj ) = -adj + j = 1,2,..., N, (10) 



fe=i 



where E nj is given by y {N+1) {T nj )P'(t nj ), r nj £ [a n _i,a n ] and P(t) = nf=o( t ~ *»»fc)- Thi s 
can be written in matrix form as 

D y-f* = -^ + jNTTfv 

where f y and E are the vectors of entries (f y )j = f(y n j,tnj) and Ej = E n j, respectively. 
From ([7j) and ffTTT) we get 



i»-{||<iid-'ii (|| / ,_ /(5 || + I J|!L Jf ) 



Since / is sufficiently smooth and satisfies a Lipschitz condition, there exist constants £ 
and /C such that 

where h = max^ \t n j — t n j-i|- This yields 

ICWD-Hh" 

\\v-£\\ < 



whenever < 1. For evenly spaced points with step h, t n j = a n _i + jh, j = 

1, 2, . . . , N, the matrix D becomes A/h, where A is a matrix whose elements are explicitly 
independent of h. Note that A is a matrix whose dimension grows as h — > 0. Therefore, 
the local error between the exact solution y n and the approximation ^ n on the nth interval 
becomes 

hn - (iv + i)(i-/ l /: n ||A-i||)- ( } 
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To get a global bound for the true error on [a, b], first let us define x n as the vector whose 
entries are the exact values of x(t) at the nodes of [a n _i,a n ], i.e., (x n )j = x(t n j). Now, 
note that f[T2"j) gives the true error \\x n — £ n || whenever the initial condition a n -x (cf. Eq. 
dSJ) for notation) is substituted by the exact value (x„_i)at. Let 5 n be the difference 

S n = (Xn)N -«n= (Xn)N - (£n)iV- 

Therefore, (11 01) becomes 

Y^D jk x nk - f(x nj ,t nj ) = -a n -i(d n )j - <y ft _i(d„)j + ^ j = 1, 2, . . . , N, 

where £ is now the Lagrange interpolation error for x(t). The same argument as above 
yields 

ilA^I! / )Ch N+1 



llXn ~ U ~ (l-fe£||A-i||) l(iVTI) + NhdM ^ ) • ^ 
where we have used the inequality ||d n || < NcIm and the definition 

d M = max \(d n )j\. 

l<j<N 
l<n<M 

By direct inspection we find that for equispaced points 

P'(t nl ) P'{t nN ) 1 



d 



M 



P'(tno)(tnl — t n o) P' (t n o) (tnN ~ t n o) Nh 



Taking into account this result and the fact that 5o = 0, the recursive use of ( fT3l) yields 
the global bound for the true error 

.. MK M \\^~ l \\h N + l 
Ft - £r < 



(iV+l)(l-^£ M ||A- 1 | 



where and £t stand for the vectors formed with the exact and approximate values at 
all the nodes t n j, j — 1, . . . , N, n — 1, . . . , M. This proves the 



Theorem 5.1. The method given by |?P is convergent and of order 0(h 

In addition, it is important to note that, as discussed by Celia and Gray pQ, a difference 
approximation as that given by the differentiation matrix (J3J) and from which SCS is 
derived, satisfies the equation 

i=o 
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imposed by the consistency requirement, if 



1 1 ■•• 1 

{to — tk) (ti — tk) ■ ■ ■ (tN — tk) 



( lk\\ 

lk2 
Tfe4 

\lkN) 



which is a straightforward consequence of the properties of the differentiation matrix (T5]). 
The analogous expression for higher order derivatives can be obtained in a similar manner. 
Therefore, the proposed method, SCS, is consistent. 

Furthermore, conditional stability of the difference scheme with respect to small changes 
on the initial value can also be easily proved. Let a and a be different initial values for a 
subinterval, and £ a , £a the corresponding solutions calculated with the method. Then 



Ma-&\\ < 



fC\\D 



-ii 



lot 



Q\ 



l-C\\D- l \ 



Therefore, small changes in the initial condition yield small changes in the numerical 
solution whenever £||Z} _1 || < 1. 



6 Numerical tests 



In this section we test SCS with some benchmark problems used frequently in the litera- 
ture. SCS is tested against two quality MATLAB solvers, ODE15s, with default param- 
eters (order 5) and quasi- const ant step size and ODE45 [TJ [S]. It should be noted that 
SCS, as used here, is not based on error control algorithms or adaptive step size strategies. 
We consider default blocks with N = 5 equispaced points in each subinterval (in order 
to compare methods of the same order) until the whole time interval is covered. In the 
following tables SCS stands for the proposed method and E is the absolute error and \\E\\ 
stands for the euclidean norm for one-dimensional problems or the Frobenius norm for 
two-dimensional problems. 

All the tests were run in MATLAB 7.2 using a personal computer with 2Gb RAM and 
Intel © processor running at 1.99GHz. 



6.1 First Order equations 
Example 1 

Let us consider the stiff problem [S] 

x' = -100x + 10, x(0) = l, (14) 
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Table 1: Results for the IVP CO} at selected points of [0,0.2]. E SC s 
and -EbDEi5s are the absolute errors obtained by using SCS and ODE15s 
respectively. The corresponding norms are H-Escsll = 0.0000714 and 
||£odei5s|| = 0.000242. 



t 


-Escs 


-EoDE15s 


0.00 


0.0 


0.0 


0.02 


0.0000688546 


0.000219504 


0.04 


0.0000186422 


0.0000679712 


0.06 


3.78549xl0- 6 


0.0000572691 


0.08 


6.83273 xl0~ 7 


0.0000134028 


0.10 


1.15621xl0~ 7 


0.0000405086 


0.12 


1.87825xl0~ 8 


4.75526 x- 6 


0.14 


2.96643xl0~ 9 


2.15198x- 7 


0.16 


4.5894xl0~ 10 


0.0000226791 


0.18 


6.9895xl0- n 


0.0000164057 


0.20 


1.0513xl0- n 


3.69695xl0~ 6 



for < x < 0.2. The solution is x(t) = (1 + 9e~ 100t )/10 and the results are shown in 
Table 1 . They are compared with those obtained with ODE15s at some points of [0, 0.2] 

Example 2 

Let us now consider the simple problem 

x = lOOx, x{0) = 1, (15) 

for < x < 0.1. Here, x(t) = e 100 *, and the results are shown in Table 2 and compared 
with those obtained with ODE15s at some points of [0,0.1] 

Example 3 

Next, we consider the nonlinear problem [9] 

x' = 5e 5 *(x - tf + 1, z(l) = -l, (16) 

for < t < 1. For this equation, x = y — e~ 5t . The results are shown in Table 3 and 
compared with those obtained with ODE45 and ODE15s at some points of [0, 1]. 
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Table 2: Results for the IVP (T5]) at selected points of [0,0.1]. 
i?SCS and -Eodei5s are the absolute errors obtained by using SCS and 
ODE15s respectively. The corresponding norms are H-Escsll = 8.03 and 
||-^ODE15s|| = 396.2. 



t 


£scs(xl0 2 ) 


-^ODElSslxlO 2 ) 


0.00 


0.0 


0.0 


0.02 


0.00000535 


0.00046360 


0.04 


0.00007917 


0.00452663 


0.06 


0.00087755 


0.04615606 


0.08 


0.00864604 


0.44160702 


0.10 


0.07986052 


3.93706426 



Table 3: Results for the IVP ([TBI) at selected points of [0,1]. Escs, 
-Eodeiss and -E0DE45 are the absolute errors obtained by using SCS, 
ODE15s and ODE45 respectively. The corresponding norms are 
liases II = 6.7 x 10- 9 , ||£odei 5s || = 6 x 10- 4 and ||£ DE45|| = 3 x 10~ 4 . 

t Escs -£<)DE15s -EODE45 

0.2 5.19952E-10 5.74683E-04 2.32239E-04 

0.4 6.99985E-11 7.58958E-05 1.74210E-04 

0.6 9.39138E-12 1.69276E-04 7.74205E-05 

0.8 1.13487E-12 6.44637E-05 3.22118E-05 

1.0 6.68797E-09 4.02615E-05 1.07195E-05 
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Table 4: Results for the IVP CED at selected points of [0,50]. E scs , 
-Eodei5s and -E0DE45 are the absolute errors obtained by using SCS, 
ODE15s and ODE45 respectively. The corresponding Frobenius norms 
are ||£ S cs|| = 1-1256 x 10~ 3 , ||E DEi5s|| = 1-7566 and ||£ode4 5 || = 
1.7558. 



t 



SCS 



E, 



ODE15s 



E, 



ODE45 



10.00 
20.00 
30.00 
40.00 
50.00 



4.35870E-04 
4.32250E-05 
2.37190E-05 
1.16350E-05 
5.35100E-06 



3.67880E-01 
1.35340E-01 
4.97870E-02 
1.83160E-02 
6.73790E-03 



3.67880E-01 
1.35340E-01 
4.97870E-02 
1.83150E-02 
6.73720E-03 



6.2 Second order problems 

In this subsection, we compare the performance of SCS against ODE15s and ODE45 in 
second order problems. 

Example 4 

Let us consider the problem (TO] 

^ = -0.1x1-199.9x2, x' 2 = -200x2, (17) 

with xi(0) = 2, £2(0) = 1, < t < 50. The solution is given by 

xi(t) = exp(-O.lt) + exp(-200t) x 2 (t) = exp(-200t), 

and the results corresponding to the first component of the solution at some points are 
shown in Table 4. 

Example 5 

Let us consider the Lotka-Volterra system 

x\ = x x (0.76 - 0.45 x 2 ), x' 2 = -x 2 (0.18 - 0.82 x x ), (18) 

with Xi(0) = 0.1, X2(0) = 0.1, < t < 1. The solver ODE45 requires 67 function 
evaluations in 11 successful steps. SCS requires 36 block function evaluations and requires 
to solve 5 linear systems. The results can be compared in Table 5. 
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Table 5: Differences |(a;i)scs - (^i)ode45| and |(x 2 )scs - (^2)ode4sI of 
the components (a?i(i), %2(i)) for the LotkaVolterra system (Tl8j) . com- 
puted by using SCS and ODE45 respectively at selected points of [0, 1]. 



t 


(^l)sCS - (xi)oDE45| 


|(>2)scs - (a:2)oDE45| 


0.25 


7.13490E-09 


1.18070E-09 


0.50 


1.68620E-08 


2.97240E-09 


0.75 


2.93810E-08 


5.61470E-09 


1.00 


4.54880E-08 


9.59720E-09 



7 Final Remarks 

The cornerstones of a simple and convergent method to solve initial value problems in 
ordinary differential equations have been presented. Due to its simplicity, it can be easily 
implemented without the need of knowing previous function values except for the initial 
value. The implementation can be made with equispaced points (as in this work) or can 
be made adaptive. 

The numerical tests show that it is indeed a competitive option, which produces accurate 
results in a wide range of problems. Only some very specific examples were discussed in 
this paper, but SCS can be easily extended to solve high order, non linear, and vector 
systems of differential equations. It is important to note that SCS is not only a compet- 
itive method, but also a innovative one, since it is a block implicit method, which makes 
it different form the standard implicit methods for the numerical solution of ordinary dif- 
ferential equations. In future papers, the application of these ideas to boundary problems 
and partial differential equations will be discussed. 
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